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(Dated: December 30, 2011) 

We study the ground state (GS) and the phase transition in a hexagonal-close-packed lattice with 
both XY and Ising models by using extensive Monte Carlo simulation. We suppose the in-plane 
interaction Ji and inter- plane interaction J2, both antiferromagnetic. The system is frustrated with 
two kinds of GS configuration below and above a critical value of = J1/J2 (t?c)- For the Ising 
case, one has rjc = 0.5 which separates in-plane ferromagnetic and antiferromagnetic states, while 
for the XY case ?7c = 1/3 separates the collinear and non coUinear spin configurations. The phase 
transition is shown to be of first (second) order for rj > (<)J7c. The spin resistivity is calculated for 
the Ising case. It shows a rounded maximum at the magnetic transition in the second-order region, 
and a discontinuity in the first-order region of r). 

PACS numbers: 05.70.Fh ; 75.10.Hk ; 75.40.Mg ;75.47.-m : 



I. INTRODUCTION 

Frustrated spin systems have been subject of intensive 
investigations during the last 30 years [ij . These systems 
are very unstable due to the competition between an- 
tagonist interactions or due to a geometrical frustration 
such as in the antiferromagnetic (AF) triangular lattice. 
A spin is said frustrated when it cannot find an orien- 
tation which " fully" satisfies all the interactions with its 
neighbors [3,11. As a consequence of the frustration, the 
ground state (GS) is very highly degenerate. In the Ising 
case the GS degeneracy is often infinite as in the triangu- 
lar lattice, the face-centered cubic (FCC) and hexagonal- 
close-packed (HCP) lattices, with AF interaction. In the 
case of vector spins, the GS is non collinear such as the 
120-degree configuration in the XY and Heisenberg AF 
stacked triangular lattice (STL). In two dimensions (2D), 
several frustrated systems with Ising spin model have 
been exactly solved [1, [H]- Among the most interest- 
ing models one can mention the frustrated generalized 
Kagome lattice Q and the honeycomb lattice where 
exotic features such as the existence of several phase tran- 
sitions, the reentrance and the disorder lines have been 
exactly found by mapping these systems into vertex mod- 
els Q. In three dimensions (3D), the situation is com- 
plicated. The renormalization group (RG) [o', 'l^, which 
provided a good understanding of the nature of the phase 
transition in non frustrated systems, encounters much of 
difficulties in application to frustrated systems. Among 
the most studied subjects during the last 20 years, one 
can mention the nature of the phase transition in the XY 
and Heisenberg STL. After a long debate [ll| on whether 
it is a second-, or a first-order transition or it belongs to a 



new universality class, the controversy has recently ended 
with the conclusion of a first-order transition 12. [l3|. 

In this paper, we are interested in the HCP antifer- 
romagnet with Ising and XY spin models. Our purpose 
is to study its properties such as the ground state, the 
phase transition and the spin transport, in the case of 
anisotropic exchange interactions. The isotropic nearest- 
neighbor (NN) AF interaction has been studied for Ising 
[1^ and XY and Heisenberg spins [l5|. These isotropic 
cases have been shown to undergo a phase transition of 
first order, and the infinite GS degeneracy is reduced to 
6 at low temperatures [13, EE] • The effect of anisotropic 
interaction and anisotropy on the GS in the case of vector 
spins has been studied [l6|. 

In section [Hi we show our model and analyze the 
ground state. Results of Monte Carlo (MC) simulations 
are shown and discussed in section IIIII Concluding re- 
marks are given in section [V] 



II. MODEL AND GROUND-STATE ANALYSIS 

We consider the HCP lattice shown in Fig. [TJ The 
stacking direction is z. The Hamiltonian is given by 



■H - - ^ Jij S j 



(1) 
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where Si is the spin at lattice site i and is the AF in- 
teraction between nearest-neighbors (NN). We suppose 
that Jij — Ji if the NN are on the xy triangular plane, 
and Jij — J2 if the NN are on two adjacent planes (see 
Fig. [T]). The GS can be determined by the steepest- 
descent method. This method is very simple [l3l (i) we 
generate an initial configuration at random (ii) we calcu- 
late the local field created at a site by its neighbors using 
([T|) (iii) we align the spin of that site along the calculated 



2 




that each tetrahedron has two up and two down 
spins. The transition between the two configura- 
tions is determined as fohows: one simply writes 
down the respective energies of a tetrahedron and 
compares them 




FIG. 1. HCP lattice. The in-plane and inter-plane interac- 
tions are indicated by Ji and J2. 



local field to minimize its energy (iv) we go to another 
site and repeat until all sites are visited: we say we make 
one sweep (v) we do a large number of sweeps per site 
until a good convergence is reached. 

However, one can also minimize the interaction energy 
as shown below to calculate the GS configuration. Since 
both interactions are AF (negative) , for simplicity, we fix 
J2 = — 1 and vary Ji. The unit of energy is taken as | J2I 
and the temperature T is in the unit of \J2\lkB where 
kB is the Boltzmann constant. 

Let us recall the GS in the case of isotropic interaction, 
namely Ji = J2 15]. For the HCP lattice, each spin is 
shared by eight tetrahedra (four in the upper half-space 
and four in the lower half-space along the z axis) and 
a NN bond is shared by two tetrahedra. The GS spin 
configuration of the system is formed by stacking neigh- 
boring tetrahedra. In the GS, one has two pairs of an- 
tiparallel spins on each tetrahedron. Their axes form an 
arbitrary angle a. The degeneracy is thus infinite (see 
Fig. 2a of Ref. [11] )■ Of course, the periodic boundary 
conditions will reduce a number of configurations, but 
the degeneracy is still infinite. Of these GS's, one partic- 
ular family of configurations of interest for both XY and 
Heisenberg cases is when a = 0. The GS is then collinear 
with two spins up and the other two down. The stack- 
ing sequence is then simplest: there are three equivalent 
configurations since there are three ways to choose the 
parallel spin pair in the original tetrahedron. 

We now examine the case where Ji 7^ J2- 

• Ising case: 

The steepest descent method with varying Ji ( J2 — 
— 1 gives two kinds of GS spin configuration: the 
first consists of xy ferromagnetic planes stacked an- 
tiferromagnetically along the z direction, while the 
second one is the stacking of xy AF planes such 



One sees that Ei < E2 when Ji > O.5J2, i.e. | Ji| < 
0.5| J2I. Thus the first configuration is more stable 
when I Ji| < 0.5| J2I. 

• XY case: 

Consider one tetrahedron of the HCP lattice. The 
Hamiltonian for this cell is given by 
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Suppose that |S,;| = 1, one has 

He — ~Ji [cos a -I- cos/3 + cos(a + /3)] 
— J2 [cos 7 -I- cos(7 — a) + cos(7 - 



(4) 



m (5) 



where the angles are defined in Fig. [51 The steep- 
est descent method shows that while /3 and a have 
unique values for a given Ji/ J2, a is arbitrary, just 
as in the case of isotropic interaction [l5| discussed 
above. To simplify the formulae, we take a = in 
the following. The energy of the cell is written as 

-He = - Ji [1 + 2cosl3] - J2 [2 cos 7 + cos(7 - /3)] (6) 

The critical values of /? and 7 are determined from 
the relations 



dp 
&Hc 
97 



= 2 Ji sin /3 - 2 J2 sin(7 - /3) = (7) 
==2J2sin7 + J2sin(7-/?) = (8) 



We find the following solutions: 

(i) /3 = 0, 7 - 0, 

(ii) /3 = 0, 7 = TT, 

(iii) /3 = TT, 7 = 0, 

(iv) /S = TT, 7 = TT, and 

(v) cos/3= 

By comparing the energy values at these solutions, 
we obtain the minimum energy with the last solu- 
tion: one has 



5 rn^^- l+3(.7i/J2)" 
4' COS7 4(Ji/j2) 
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1 + 3(Ji/ J2)^ 
2(Jl/J2) 



(9) 
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FIG. 2. Ground state in the XY case. The tetrahedron is 
projected on the xy plane. The spins are numbered from 1 to 
4. See text for comments. 
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FIG. 3. Ground state in the XY case (color online). The 
angles /? (black circles) and 7 (blue void circles) and their 
cosinus are shown as functions of -q — J1/J2. Non coUinear 
GS configurations occur in the region 1/3 < r; < 1 . See text 
for comments. 



where cos (3 and cos 7 are given above. 

Because — 1 < cos/3 < 1 and — 1 < cos 7 < 1, the 
above solution is valid for | < J1/J2 < 1. We plot 
cos^, cos 7, /3 and 7 in Fig. [3] where we observe 
that the non coUinear GS configuration occurs in 
the interval 1/3 < i] = Ji/ J2 < 1. 



interaction IJ2I = 1 is used as the unit of energy. We 
use here MC simulations with a histogram technique to 
detect first-order transition. We equilibrate the system 
and average physical quantities over several millions of 
MC steps per spin. The averaged energy (U) and the 
heat capacity Cy are calculated by 



([/) = in) 



(10) 
(11) 



where (...) indicates the thermal average taken over mi- 
croscopic states at T. 

The order parameter M is defined from the sublattice 
magnetization by 



(12) 



K ieK 



where belongs to the sublattice K. Note that there are 
at least four sublattices which define the ordering of the 
spins on the tetrahedra if a long-range order is observed. 
The susceptibility is defined by 



X 



(13) 



In MC simulations, we work at finite sizes, so for each 
size we have to determine the " pseudo" transition which 
corresponds in general to the maximum of the specific 
heat or of the susceptibility. The maxima of these quan- 
tities need not to be at the same temperature. Only 
at the infinite size, they should coincide. The theory 
of finite-size scaling fl8|-[20| permits to deduce properties 
of a system at its thermodynamic limit. In order to 
check the first-order nature of the transition, we used 
the histogram technique which is very efficient in detect- 
ing weak first-order transitions and in calculating critical 
exponents of second-order transitions. [Tol [20| The main 
idea of this technique is to make an energy histogram at 
a temperature Tq as close as possible to the transition 
temperature. Often, one has to try at several tempera- 
tures in the transition region. Using this histogram in the 
formulae of statistical physics for canonical distribution, 
one obtains energy histograms in a range of temperature 
around Tq. In second-order transitions, these histograms 
are gaussian. They allows us to calculate averages of 
physical quantities as well as critical exponents using the 
finite-size scaling. In first-order transitions, the energy 
histogram shows a double-peak structure. 



III. PHASE TRANSITION: RESULTS 

We consider a sample size oi Lx Lx where L and 
vary from 12 to 36 but can be different from L in order 
to detect the dependence of the GS on Ji . The exchange 



A. Ising case 

In the following, the results in different figures are 
shown with L = 18 and = 8 (16 atomic planes along 
z). Since the GS changes at tJc = 0.5, we show here ex- 
amples on both sides of this value. Figured shows the en- 
ergy per spin E, the specific heat per spin Cy, the order 
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FIG. 4. Ising case: Energy specific heat Cv, order param- 
eter M and susceptibility x versus T for r\ = J1/J2 = 0.3. 
See text for comments. 



FIG. 5. Ising case (color online): Energy E and order param- 
eter M for r) = J1/J2 = 0.85 (blue void circles) and 1 (red 
triangles). See text for comments. 



P(E) 




FIG. 6. Ising case (color online): Energy histogram P{E) 
versus E for t; = 1( red triangles), 0.85 (blue void circles) and 
0.3 (black circles). See text for comments. 



B. XY case 



parameter M and the susceptibility x, for 77 = 0.3. The 
transition is of second order. On the other side, we show 
in Fig. [5] the energy per spin and the order parameter 
versus T, for 77 = 0.85 and 1. We find a strong first-order 
transition in both cases. The discontinuity of E and M 
at the transition is very large. We show in Fig. [Slthe en- 
ergy histogram taken at the transition temperature for 
three values ry = 0.3 (black), 0.85 (blue) and 1 (red). As 
seen, the first case is a gaussian distribution indicating a 
second-order transition, while the last two cases show a 
double-peak structure indicating a first-order transition. 

We have calculated the critical temperature Tc as a 
function of 77. The phase diagram is shown in Fig. [7] 
where I and II indicate the ordering of the first, and sec- 
ond kinds, respectively. P indicates the paramagnetic 
phase. Note that the transition line between I and P 
is a second-order line, while that between II and P is a 
first-order line. 



In the XY case, the change of the GS takes place at 
r] — 1/3. Let us show in Fig. [8] the result for 77 = 0.3 
where the GS is composed of ferromagnetic planes anti- 
ferromagnetically stacked in the z direction. The transi- 
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FIG. 7. Ising case: Tc versus 77. I, II and P denote the first, 
second and paramagnetic phases, respectively. See text for 
comments. 
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FIG. 8. XY case: Energy E, specific heat Cv, order parame- 
ter M and susceptibility x versus T for i] — 0.3. See text for 
comments. 



tion is of second order. 

We show now the result for the non colUnear GS region 
in Fig. |9l The energy and the order parameter show 
clearly a discontinuity at the transition for i] = 0.58 and 
1. Using the histogram method described above, we have 
calculated the histogram shown in Fig. [10] for 77 = 0.3, 
0.58 and 1. For ry = 0.3 which is in the coUinear region of 
the GS, the histogram is gaussian, confirming the second- 
order transition observed in the data shown above. For 
T] — 0.58 and 1 belonging to the non coUinear region, the 
histogram shows a two-peak structure which confirms the 
first-order character of the transitions in this region. 

The two peaks are very well separated with the dip 
going down to zero, indicating an energy discontinuity. 
The distance between the two peaks is the latent heat 
AE. 

We show in Fig. [11] the transition temperature versus 
r] where I and II indicate the coUinear and non coUinear 
phases, respectively. P denotes the paramagnetic state. 




FIG. 9. XY case (color online): E and A/ versus T for ri = 
0.58 (blue void circles) and 1 (red triangles). See text for 
comments. 




FIG. 10. XY case (color online): Energy histogram P versus 
E for Tj — 0.3 (black circles), 0.58 (blue void circles), 1 (red 
triangles) at the respective transition temperatures. See text 
for comments. 



The line separating I and P is a second-order transition 
line, while that separating II and P is the first-order one. 

To close this section, we emphasize that all 3D frus- 
trated systems known so far undergo a first-order tran- 
sition: let us mention the AF STL [H, [11], the FCC 
antiferromagnets pTj . the simple cubic fully frustrated 
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FIG. 11. XY case: To versus 77. I, II and P denote the 
coUinear, non coUinear and paramagnetic phases, respectively. 
See text for comments. 
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lattices [22-1251, helimagnets [261 and the HCP lattice 
studied here. 



IV. SPIN TRANSPORT IN THE ISING CASE 

We have recently studied the spin transport in mag- 
netically ordered materials by using an efficient MC sim- 
ulation method [13, . The tour-de- force of the method 
was realized when we studied the semiconducting MnTe 
where the agreement with experiment is excellent [29] . In 
ferromagnets, the spin resistivity shows a sharp peak at 
the magnetic phase transition. The origin of this anomaly 
comes from the spin-spin correlation |30l432l |. In antifer- 
romagnets, such a peak is not sharp and often a broad 
maximum [33[. We have also studied the spin resistivity 
in two frustrated cases: the FCC Ising case [3J| and the 
Ji — J2 simple cubic lattice [35|. The results show that 



the discontinuous transition in these systems gives rise 
to a discontinuity of the spin resistivity. 

Using the method which has been described in detail 
elsewhere (27,, ^SJ,] , we calculate the current of itinerant 
spins moving across the crystal under an electric field e 
applied in the x direction. We suppose that each itinerant 
spin CTi interacts with the lattice spins around it within 
a sphere of radius Di at any point on its trajectory. The 
interaction Hamiltonian reads 



(14) 



where the sum is performed over all lattice spins in the 
sphere around the itinerant spin ai and lij is the in- 
teraction which depends on the distance between and 
Sj . We suppose also the following interaction between an 
itinerant spin with its neighboring itinerant spins within 
a sphere of radius D2 



Hi — — KijCJi 



For simplicity, we suppose 



hj = /oexp(-Bry) 
Kij = KoeyiY>{-Crij) 



(15) 



(16) 
(17) 



where /q, B, Kq and C are constants chosen in such a 
way that the energy of an itinerant spin is smaller than 
that of a lattice spin. This choice is made to avoid the 
influence of itinerant spins on the ordering of the lat- 
tice spins. Note that in the almost-free electron model, 
Iq — Kq — 0, and in semiconductors they are larger but 
still weak with respect to Ji and J2 . A discussion on the 
choice of different parameters has been given for example 
in Ref. [27| . We suppose here a concentration of one itin- 
erant per two lattice cells. At such a low concentration, 
the averaged distance between itinerant spins is much 
larger than the cutoff distance D2- In spite of this, due 



to the attractive nature of their interaction, Eq. (IT7|) . we 
have to use a chemical potential term to insure that itin- 
erant spins do not form clusters. This potential is taken 
of the form D[n(r) — uq] where D is a positive constant, 
n(f) the concentration of itinerant spins in the sphere of 
cutoff radius D2 centered at the position r of the itinerant 
spin under consideration, and uq the averaged concentra- 
tion. We perform the simulation by takin g in to account 
the relaxation time of the lattice spins l28l [29| which 
is given by 



A 



TL 



\l-T/Tc 



(18) 



where y4 is a constant, v the correlation critical exponent, 
and z the dynamic exponent. We see that as T tends 
to Tc, tl diverges. This phenomenon is known as the 
critical slowing-down. For the Ising spin model, v — 
0.638 (3D Ising universality) and z = 2.02 [s^. We have 
previo usly shown that tl strongly affects the shape of p 
at Tc nil. By choosing A = 1, we fix ti, = 1 at T = 2Tc 
deep inside the paramagnetic phase far above Tc- This 
value is what we expect for thermal fluctuations in the 
disordered phase. 

The spin resistivity p versus T for the two typical cases 
rj — 0.3 and 1 is shown in Fig. [T^] where distances Di 
and D2 are in unit of the NN distance, energy constants 
/o, Kq and D are in unit of IJ2I = 1. We observe here 
that in the second-order region p has a rounded maxi- 
mum and in the first-order region it undergoes an almost 
discontinuous jump at the transition. 

We have varied the radius Di to see its effect on the 
resistivity p at the transition. We observe the same effect 
which has been seen in some other antiferromagnets [l^, 
I35I : at a given T, p oscillates slightly with distance. We 
have found that this oscillation comes from the oscillatory 
behavior of the difference between the numbers of up and 
down spins in the sphere as Z^i varies. 



V. CONCLUSION 

We have studied in this paper some properties of the 
HCP antiferromagnet with Ising and XY spin models. 
The in-plane Ji and inter-plane J2 interactions are sup- 
posed to be different. As a result, the GS spin configu- 
ration depends on the ratio 77 = Ji/J2- We show that 
there exists a critical value rjc where the GS changes. 
For the Ising case, we find r^c = 0.5 below (above) which 
the spins in the xy planes are ferromagnetic (antiferro- 
magnetic). For the XY case, the GS is coUinear below 
?7c = 1/3, and is non coUinear above that value. The 
nature of the transition changes from a second order be- 
low rjc to a first order above rjc for both Ising and XY 
cases. We have also studied the spin resistivity in the 
Ising case. We found that the shape of p depends on the 
nature of the transition: in the second-order region, a 
rounded maximum is observed at the transition while in 
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FIG. 12. Ising case: Spin resistivity p versus T for rj = 0.3 
(upper) and 1 (lower). L — 18, Lz — 8 (16 planes in the z 
direction), Di = Da = 1, e = 1, /o = 2, Kq = 0.5, B = 
C = 1, A = 1, D = 0.5. All distances are in unit of the NN 
distance, energy constants are in unit of jJaj = 1. See text 
for comments. 



the first-order region, p undergoes a discontinuity at the 
transition as we have observed in other frustrated cases. 
These findings may help to understand the transition na- 
ture and the spin transport in different compounds with 
HCP structure. 
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